rm(list = ls())
require(miceadds)
require(stats4)
require(lmtest)
require(bbmle)
require(msm)
require(nlWaldTest)
require(sandwich)
require(stargazer)
require(mfx)
require("multiwayvcov")
require(car)
require(alr3)
 
#data2old=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/infoexp2_data/responses.csv", header=TRUE)
#data13=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_july25_type13clean.csv", header=TRUE)
#data15=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_August12_type15(clean).csv", header=TRUE)

data2old=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/infoexp2_data/responses.csv", header=TRUE)
data13=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_july25_type13clean.csv", header=TRUE)
data15=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_August12_type15(clean).csv", header=TRUE)


dataframe2old=data.frame(data2old)
dataframe13=data.frame(data13)
dataframe15=data.frame(data15)


#fancy regression for easy clustering by Markus Conrad https://www.r-bloggers.com/2021/05/clustered-standard-errors-with-r/
ols <- function(form, data, robust=FALSE, cluster=NULL,digits=3){
  r1 <- lm(form, data)
  if(length(cluster)!=0){
    data <- na.omit(data[,c(colnames(r1$model),cluster)])
    r1 <- lm(form, data)
  }
  X <- model.matrix(r1)
  n <- dim(X)[1]
  k <- dim(X)[2]
  if(robust==FALSE & length(cluster)==0){
    se <- sqrt(diag(solve(crossprod(X)) * as.numeric(crossprod(resid(r1))/(n-k))))
    res <- cbind(coef(r1),se)
  }
  if(robust==TRUE){
    u <- matrix(resid(r1))
    meat1 <- t(X) %*% diag(diag(crossprod(t(u)))) %*% X
    dfc <- n/(n-k)    
    se <- sqrt(dfc*diag(solve(crossprod(X)) %*% meat1 %*% solve(crossprod(X))))
    vcov <- solve(crossprod(X)) %*% meat1 %*% solve(crossprod(X))
    res <- cbind(coef(r1),se)
    }
  if(length(cluster)!=0){
    clus <- cbind(X,data[,cluster],resid(r1))
    colnames(clus)[(dim(clus)[2]-1):dim(clus)[2]] <- c(cluster,"resid")
    m <- dim(table(clus[,cluster]))
    dfc <- (m/(m-1))*((n-1)/(n-k))
    uclust  <- apply(resid(r1)*X,2, function(x) tapply(x, clus[,cluster], sum))
    se <- sqrt(diag(solve(crossprod(X)) %*% (t(uclust) %*% uclust) %*% solve(crossprod(X)))*dfc)   
    vcov <- solve(crossprod(X)) %*% (t(uclust) %*% uclust) %*% solve(crossprod(X))
    res <- cbind(coef(r1),se)
    }
  res <- cbind(res,res[,1]/res[,2],(1-pnorm(abs(res[,1]/res[,2])))*2)
  res1 <- matrix(as.numeric(sprintf(paste("%.",paste(digits,"f",sep=""),sep=""),res)),nrow=dim(res)[1])
  rownames(res1) <- rownames(res)
  colnames(res1) <- c("Estimate","Std. Error","t value","Pr(>|t|)")
  out=list(res1,vcov)
  return(out)
}

priors=c(0.5,0.6,0.75,0.85)

expframe=rbind(dataframe13,dataframe15,dataframe2old)
expframe=subset(expframe,uid>200)

#note that player 234 did not finish and so has no real data
expframe=subset(expframe,uid!=234)


#definitions
expframe$priorA=(expframe$qid==1)*0.5+(expframe$qid==5)*0.6+(expframe$qid==6)*0.75+(expframe$qid==7)*0.85
expframe$priorB=1-expframe$priorA
expframe$Aresponse=(expframe$response==6)
expframe$Bresponse=(expframe$response==7)
expframe$stateA=(expframe$state==4)
expframe$stateB=(expframe$state==10)
expframe$correct=expframe$stateA*expframe$Aresponse+expframe$stateB*expframe$Bresponse
uidlist=unique(expframe$uid)

#find N
#**
length(uidlist)
#**

#Initialize vectors
pAgiven2=vector('numeric')
restrictionpAgiven1=vector('numeric')
pAgiven1=vector('numeric')
prob=vector('numeric')
sterr=vector('numeric')
tstat=vector('numeric')

#check for violations
NIASindiframe=data.frame(uid=vector('numeric'),pval0.5=vector('numeric'),pval0.6=vector('numeric'),pval0.75=vector('numeric'),pval0.85=vector('numeric'))
for(i in 1:length(uidlist)){
playframe=subset(expframe,uid==uidlist[i])
frame50=subset(playframe,priorA==0.5)
output50=ols(Aresponse~stateA,data=frame50,robust=TRUE, cluster=NULL)
model50=output50[[1]]
vcov50=output50[[2]]
pAgiven2[1]=model50[1]
restrictionpAgiven1[1]=((2*priors[1]-1)/priors[1])+pAgiven2[1]*((1-priors[1])/priors[1])
pAgiven1[1]=model50[1]+model50[2]
sterr=as.numeric(deltamethod(~x1+x2-((2*0.5-1)/0.5)+x1*((1-0.5)/0.5),model50[1:2],vcov50))
tstat=(pAgiven1[1]-restrictionpAgiven1[1])/sterr
prob[1]=pnorm(tstat)


frame60=subset(playframe,priorA==0.6)
output60=ols(Aresponse~stateA,data=frame60,robust=TRUE)
model60=output60[[1]]
vcov60=output60[[2]]
pAgiven2[2]=model60[1]
restrictionpAgiven1[2]=((2*priors[2]-1)/priors[2])+pAgiven2[2]*((1-priors[2])/priors[2])
pAgiven1[2]=model60[1]+model60[2]
sterr=as.numeric(deltamethod(~x1+x2-((2*0.6-1)/0.5)+x1*((1-0.6)/0.5),model60[1:2],vcov60))
tstat=(pAgiven1[2]-restrictionpAgiven1[2])/sterr
prob[2]=pnorm(tstat)



frame75=subset(playframe,priorA==0.75)
output75=ols(Aresponse~stateA,data=frame75,robust=TRUE)
model75=output75[[1]]
vcov75=output75[[2]]
pAgiven2[3]=model75[1]
restrictionpAgiven1[3]=((2*priors[3]-1)/priors[3])+pAgiven2[3]*((1-priors[3])/priors[3])
pAgiven1[3]=model75[1]+model75[2]
sterr=as.numeric(deltamethod(~x1+x2-((2*0.75-1)/0.5)+x1*((1-0.75)/0.5),model75[1:2],vcov75))
tstat=(pAgiven1[3]-restrictionpAgiven1[3])/sterr
prob[3]=pnorm(tstat)


frame85=subset(playframe,priorA==0.85)
output85=ols(Aresponse~stateA,data=frame85,robust=TRUE)
model85=output85[[1]]
vcov85=output85[[2]]
pAgiven2[4]=model85[1]
restrictionpAgiven1[4]=((2*priors[4]-1)/priors[4])+pAgiven2[4]*((1-priors[4])/priors[4])
pAgiven1[4]=model85[1]+model85[2]
sterr=as.numeric(deltamethod(~x1+x2-((2*0.85-1)/0.5)+x1*((1-0.85)/0.5),model85[1:2],vcov85))
tstat=(pAgiven1[4]-restrictionpAgiven1[4])/sterr
prob[4]=pnorm(tstat)

addframe=data.frame(uid=uidlist[i],pval0.5=prob[1],pval0.6=prob[2],pval0.75=prob[3],pval0.85=prob[4])
NIASindiframe=rbind(NIASindiframe,addframe)
}

#need to clean an NaN from first part
cleanpval0.5= NIASindiframe$pval0.5[!is.na(NIASindiframe$pval0.5)]

unexpectedchangers0.5=100*sum(cleanpval0.5<0.05)/length(NIASindiframe$uid)
unexpectedchangers0.6=100*sum(NIASindiframe$pval0.6<0.05)/length(NIASindiframe$uid)
unexpectedchangers0.75=100*sum(NIASindiframe$pval0.75<0.05)/length(NIASindiframe$uid)
unexpectedchangers0.85=100*sum(NIASindiframe$pval0.85<0.05)/length(NIASindiframe$uid)

violateframe=data.frame(prior=priors,persigvio=c(unexpectedchangers0.5,unexpectedchangers0.6,unexpectedchangers0.75,unexpectedchangers0.85))
violateframe=t(violateframe)

#***********************************
stargazer(violateframe,summary=FALSE,colnames=FALSE)
#************************************

